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Abstract. 

In this article, we review some of our approaches to granular dynamics, now well 
known £Q to consist of both fast and slow relaxational processes. In the first case, 
grains typically compete with each other, while in the second, they cooperate. A typical 
result of cooperation is the formation of stable bridges, signatures of spatiotemporal 
inhomogeneities; we review their geometrical characteristics and compare theoretical 
results with those of independent simulations. Cooperative excitations due to local 
density fluctuations are also responsible for relaxation at the angle of repose; the 
competition between these fluctuations and external driving forces, can, on the other 
hand, result in a (rare) collapse of the sandpile to the horizontal. Both these features 
are present in a theory reviewed here. An arena where the effects of cooperation versus 
competition are felt most keenly is granular compaction; we review here a random 
graph model, where three-spin interactions are used to model compaction under 
tapping. The compaction curve shows distinct regions where 'fast' and 'slow' dynamics 
apply, separated by what we have called the single-particle relaxation threshold. In the 
final section of this paper, we explore the effect of shape - jagged vs. regular - on 
the compaction of packings near their jamming limit. One of our major results is 
an entropic landscape that, while microscopically rough, manifests Edwards' flatness 
at a macroscopic level. Another major result is that of surface intermittency under 
low-intensity shaking. 
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1. INTRODUCTION 

Matter in the jammed state has become a focus of interest for physicists in recent 
years. Two prime examples of this in the context of natural systems are glasses E] 
and densely packed granular media 0; while the mechanisms of jamming in each case 
show strong similarities, the ineffectiveness of temperature as a dynamical motor in 
granular media leads to vastly more surprising effects. A direct consequence of such 
athermal behaviour in sandpiles is the stable formation of cooperative structures such 
as bridges [Sj, or indeed the very existence jH] of an angle of repose [7j; neither would 
be possible in the presence of Brownian motion. We discuss our studies of these two 
effects in the first two sections of this article. The last two sections, with their focus 
on jamming, unify aspects of glasses and granular media; one of them uses random 
graphs to illustrate competitive and cooperative effects in granular compaction 8J. The 
other concerns itself mainly with the fast dynamics in the boundary layer of a granular 
column, as the jamming limit is approached; it makes clear how asymmetric grains can 
orient themselves suitably so as to waste less space, when compelled so to do 

2. On bridges in sandpiles - an overarching scenario 

The athermal nature of granular media results in the following fact: all granular 
dynamics is the result of external stimuli. These result in grains competing with each 
other to fall under gravity to a point of stability; when instead, the process is one of 
cooperation so that two or more grains fall together to rest on the substrate with mutual 
support, bridges [S] are formed. These can be stable for arbitrarily long times, since 
the Brownian motion that would dissolve them away in a liquid is absent in sandpiles 
- grains are simply too large for the ambient temperature to have any effect. As a 
result, bridges can affect the ensuing dynamics of the sandpile; a major mechanism 
of compaction is the gradual collapse of long-lived bridges in weakly vibrated granular 
media, resulting in the disappearance of the voids that were earlier enclosed JU] • Bridges 
are also responsible for jamming in granular processes, for example, as grains flow out 
of a hopper [Oj. 

We first define a bridge in more quantitative terms. Consider a stable packing of 
hard spheres under gravity, in three dimensions. Each particle typically rests on three 
others which stabilise it, in the sense that downward motion is impeded. A bridge 
is a configuration of particles in which the three-point stability conditions of two or 
more particles are linked; that is, two or more particles are mutually stabilised. Bridges 
thus cannot be formed sequentially, but are ubiquitous in generic powders. While it 
is impossible to determine bridge distributions uniquely from a distribution of particle 
positions, we are able via our algorithm to obtain the most likely positions of bridges 
in a given scenario [Sj. 

We now distinguish between linear and complex bridges via a comparison of 
Figures C] arid |2j Figure illustrates a complex bridge, i.e., a mutually stabilised cluster 



Competition and cooperation 



3 



of five particles (shown in green), where the stability is provided by six stable base 
particles (shown in blue). Of course the whole is embedded in a stable network of grains 
within the sandpile. Also shown is the network of contacts for the particles in the bridge: 
we see clearly that three of the particles each have two mutual stabilisations. Figure El 
illustrates a seven particle linear bridge with nine base particles. This is an example 
of a linear bridge. The contact network shows that this bridge has a simpler topology 
than that in Figure ^ Here, all of the mutually stabilised particles are in sequence, as 
in a string. A linear bridge made of n particles therefore always rests on = n + 2 
base particles. For a complex bridge of size n, the number of base particles is reduced 
(rib < n + 2), because of the presence of loops in their contact networks. 




Figure 1. A five particle complex bridge, with six base particles (left), and the 
corresponding contact network (right). Thus n = 5 and n^ = 6 < 5 + 2. 

An important point to note is that bridges can only be formed sustainably in 
the presence of friction; the mutual stabilisations needed would be unstable otherwise! 
Although our Monte Carlo simulations (described below) do not contain friction 
explicitly, our configurations indirectly include this: in particular, the coordination 
numbers lie in a range consistent with the presence of friction [TlH ITU IT2"]. 

2.1. Simulation details 

We have examined bridge structures in hard assemblies that are generated by a 
non-sequential restructuring algorithm ^U|, whose main modelling ingredients involve 
stochastic grain displacements and collective relaxation from them. 

This algorithm restructures a stable hard sphere deposit in three distinct stages. 

(1) The granular assembly is dilated in a vertical direction (with free volume being 
introduced homogeneously throughout the system), and each particle is given 
a random horizontal displacement; this models the dilation phase of a vibrated 
granular medium. 
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Figure 2. A seven particle linear bridge with nine base particles (left), and the 
corresponding contact network (right). Thus n = 7 and n^ = 9 = 7 + 2. 

(2) The assembly is compressed in a uniaxial external field representing gravity, using 
a low-temperature Monte Carlo process. 

(3) Individual spheres in the assembly are stabilised using a steepest descent 'drop and 
roll' dynamics to find local potential energy minima. 

Steps (2) and (3) model the quench phase of the vibration, where particles relax to 
locally stable positions in the presence of gravity. Crucially, during the third phase, the 
spheres are able to roll in contact with others; mutual stabilisations are thus allowed 
to arise, mimicking collective effects. The final configuration has a well-defined contact 
network where each sphere is supported by a uniquely defined set of three other spheres. 

The simulation method recalled above builds a sequence of static packings. Each 
new packing is built from its predecessor by a random process and the sequence achieves 
a steady state, where structural descriptors such as the mean packing fraction and 
the mean coordination number fluctuate about well-defined mean values. The steady- 
state mean volume fraction $ typically evolves to values in the range $ ~ 0.55 - 0.61, 
depending on the shaking amplitude; the mean coordination number is always Z « 
4.6 ± 0.1. We recall that for frictionless (isostatic) packings in d dimensions, Z — 2d 
(= 6 for d — 3) while for frictional packings the minimal coordination number is 
Z = d + 1 (= 4 for d = 3) [TT] : our configurations thus clearly correspond to those 
generated in the presence of friction. This is confirmed by the results of molecular 
dynamics simulations of sphere packings in the limit of high friction, which yield a 
mean coordination number slightly above 4.5 |12j . 

Each of our configurations includes N tot ~ 2200 particles. Segregation is avoided 
by choosing monodisperse particles: a rough base prevents ordering. A large number of 
restructuring cycles is needed to reach the steady state for a given shaking amplitude: 
about 100 stable configurations (picked every 100 cycles in order to avoid correlation 
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effects) are analysed, corresponding to $ = 0.56 and $ = 0.58. From these 
configurations, and following specific prescriptions, our algorithm identifies bridges as 
clusters of mutually stabilised particles 0. 

Figure H3 illustrates two characteristic descriptors of bridges used in this work. The 
main axis of a bridge is defined using triangulation of its base particles as follows: 
triangles are constructed by choosing all possible connected triplets of base particles, 
and the vector sum of their normals is defined to be the direction of the main axis of 
the bridge. The orientation angle G is defined as the angle between the main axis and 
the z-axis. The base extension b is defined as the radius of gyration of the base particles 
about the z-axis; note that this is distinct from the radius of gyration about the main 
axis of the bridge. 




Figure 3. Definition of the angle and the base extension b of a bridge. 
The main axis makes an angle with the z-axis; the base extension b is the 
projection of the radius of gyration of the bridge on the x-y plane. 



2.2. Bridge sizes and diameters: when does a bridge span a hole? 

In the following, we present statistics for both linear and complex bridges. While we 
recognise that bridge formation is a collective dynamical process, we adopt an ergodic 
viewpoint ^3] here. Inspired by polymer theory ^Hj, we visualise a linear bridge as 
a random chain which grows as a continuous curve, i.e. 'sequentially' in terms of its 
arc length s. (For complex bridges, this simplification is not possible in general - a 
direct consequence of their branched structure). This replacement of what is in reality 
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a collective phenomenon in time by a random walk in space is somewhat analogous to 
the 'tube model' of linear polymers [T^]: both are simple but efficient effective pictures 
of very complex problems. 

We first address the question of the length distribution of linear bridges. We define 
the length distribution f n as the probability that a linear bridge consists of exactly n 
spheres. We make the simplest and the most natural assumption that a bridge of size 
n remains linear with some probability p < 1 if an (n + l) th sphere is 'added' to it: this 
leads to the exponential distribution 

fn = (l-p) P n . (2.1) 

The exponential distribution above can also be derived by means of a continuum 
approach. Here, a linear bridge is viewed as a continuous random curve or 'string', 
parametrised by the arc length s from one of its endpoints. We assume also that such 
a bridge disappears at a constant rate a per unit length, either by changing from linear 
to complex or by collapsing. The survival probability S(s) of a linear bridge upto 
length s thus obeys the rate equation S = —aS and falls off exponentially, according to 
S(s) = exp(— as). Consequently, the probability distribution of the length s of linear 
bridges reads f(s) = —S(s) = a exp(— as), a continuum analogue of (|2.1j) . 

This is in good accord with the results of independent simulations, which exhibit 
an exponential decay of linear bridges of the form (j2.1j) . with a ~ 0.99 [3], which is 
clearly seen until n ~ 12. Around n ~ 8, complex bridges begin to predominate; these 
have size distributions which show a power-law decay: 

fn ~ n~ T . (2.2) 

with r « 2 0. 

We have also measured the diameter R n of linear and complex bridges of size 
n, which is such that R\ is the mean squared end-to-end distance. Our data on 
diameters and size distributions j^j indicate that linear bridges in three dimensions start 
off as planar self-avoiding walks, which eventually collapse onto each other because of 
vibrational effects; on the other hand, complex bridges look like 3d percolation clusters. 

Another issue of interest to us is the jamming potential of a bridge. A measure 
of this, in the case of a linear bridge, is its the base extension b (see Figure EI); this is 
the horizontal projection of the 'span' of the bridge. Our simulation results i5j indicate 
that three-dimensional bridges of a given length have a fairly characteristic horizontal 
extension, making it relatively easy to predict whether or not they would 'jam' a given 
hole. 

In order to compare our simulations with experiment JHlEj, we plot in Figure|3]the 
logarithm of the probability distribution of base extensions p(b) against the (normalised) 
base extension b/(b). This figure emphasises the exponential tail of the distribution 
function, and also shows that bridges with small base extensions are unfavoured. We 
note that this long tail is characteristic of three-dimensional experiments on force chains 
in granular media ^7j. The sharp drop at the origin as well as the long tail in Figure 0J 
are observed in normal force distributions obtained via molecular dynamics simulations 
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of particle packings ^H] , i n the limit of strong deformations. Realising that the measured 
forces propagate through chains of particles, we use this similarity to suggest that bridges 
are really just long-lived force chains, which have survived despite strong deformations. 
We suggest also that with the current availability of 3D visualisation techniques such 
as NMR PHj, bridge configurations might be an easily measurable and effective tool to 
probe inhomogeneities in shaken sand. 




0.0 0.5 



1.0 1.5 

b/<b> 



2.0 2.5 



Figure 4. Distribution of base extensions of bridges, for <]? = 0.58. The 
logarithm of the normalised probability distribution is plotted as a function 
of the normalised variable b/(b), where (b) is the mean extension of bridge 
bases. 



2. 3. Turning over at the top; how linear bridges form domes 

Recall that a linear bridge is modelled as a continuous curve, parametrised by its arc 
length s. We here focus on its most important degree of freedom, the tilt with respect 
to the horizontal; the azimuthal degree of freedom is neglected. Accordingly, we define 
the local or link angle 9{s) between the direction of the tangent to the bridge at point s 
and the horizontal, and the mean angle made by the bridge from its origin up to point s, 
also with the horizontal: 

e(s) = - [ S 0(u)du. (2.3) 
s Jo 

The local angle 9(s) so defined may be either positive or negative; it can even change 
sign along the random curve which represents a linear bridge. Of course, the orientation 
angle measured in our numerical simulations is positive by construction, being defined 
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as the angle between the main bridge axis and the z-axis (see Figure |3J). (Note that by 
simple geometry, this 'zenith angle' made by the bridge axis with the vertical, equals 
the mean angle made by the basal plane of the bridge with the horizontal) . 

Our simulations show that the mean angle 0(s) typically becomes smaller and 
smaller as the length s of the bridge increases. Small linear bridges are almost never 
flat as they get longer, assuming that they still stay linear, they get 'weighed down', 
arching over as at the mouth of a hopper [Hj. Thus, in addition to our earlier claim that 
long linear bridges are rare, we claim further here that (if and) when they exist, they 
typically have flat bases, becoming 'domes'. 

We use these insights to write down equations to investigate the angular distribution 
of linear bridges. These couple the evolution of the local angle 6(s) with local density 
fluctuations 0(s) at point s (with ' denoting a derivative with respect to s): 

e' = - a e-b(j) 2 + A irh (s), (2.4) 

0' = -c0 + A 2 77 2 (s). (2.5) 

The effects of vibration on each of 9 and are represented by two independent white 
noises r/i(s), r] 2 (s), such that 

faOOTfoCs')) = 2 SijSis-s'), (2.6) 

whereas the parameters a, A 2 are assumed to be constant. 

The phenomenology behind the above equations is the following: the evolution 
of 8(s) is caused, in our effective picture, by the sequential addition of particles to 
the bridge at its ends. The fluctuations of local density at a point s are caused by 
collective particle motion j20j. The first terms on the right-hand side of (12 .4j) . ()2.5j) say 
that neither 9 nor is allowed to be arbitrarily large. Their coupling via the second 
term in ()2.4|) arises as follows: if there are density fluctuations 2 of large magnitude 
at the tip of a bridge, these will, to a first approximation, 'weigh the bridge down', i.e., 
decrease the angle 6 locally. 

Reasoning as above, we therefore anticipate that for low-intensity vibrations and 
stable bridges, both density fluctuations 0(s) and link angles 9(s) will be small. 
Accordingly, we linearise ()2.4|) . obtaining thus an Ornstein-Uhlenbeck equation: 

0> = - a 6 + Air7i(s). (2.7) 

Let us make the additional assumption that the initial angle 8q, i- e -> that observed for 
very small bridges, is itself Gaussian with variance <Jq = (6q). The angle 8(s) is then a 
Gaussian process with zero mean for any value of the length s. Its correlation function 
can be easily evaluated to be |2*T] : 

(O(s)0(s')} = <e-*l-'l + {o-l - <)e-^'). (2.8) 

It follows from this that the variance of the link angle is: 

<0 2 >(*) = < + K 2 -<)e- 2as , (2.9) 

We see from the above that orientation correlations decay with a characteristic length 
given by £ = 1/a; also, in the limit of an infinite bridge, the variance (6 2 ) relaxes to 
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a 2 = — - ED- Thus, as the chain gets longer, the variance of the link angle relaxes from 
q a 

its initial value of a 2 (i.e. that for the initial link) to a 2 q for infinitely long chains. 

Given the above, it can be shown that the mean angle 0(s) will also have a Gaussian 
distribution. Its variance can be derived by inserting (j2.8J) into ()2.3j) : 

<0 2 >(s) = 2< ag-1 2 + 2 e " aS + {al - <) (1 Xr )2 - ( 2 -10) 



'"'/^■i ■ -",,1 

as as 

The asymptotic result 



(6 2 )( s )-^a-^i, (2.11) 
as as 

confirms our earlier statement that the longest bridges form domes, i.e. they have bases 
that are almost flat. Each such bridge can be viewed as consisting of a large number 
as = s/£ ^> 1 of independent 'blobs' of length £; this result suggests, yet again, strong 
analogies between linear bridges and linear polymers [TK] . 

The result (12.11}) has another interpretation. As @(s) is small with high probability 
for a very long bridge, its extension in the vertical direction reads approximately 

Z = z(s) - z(0) « s9(s), (2.12) 

so that (Z 2 ) « s 2 (Q 2 )(s) ~ 2(Ax/a) 2 s. Switching back to a discrete picture of an n-link 
chain, we have 

Z n ~ n 1/2 . (2.13) 

The vertical extension of a linear bridge is thus found to grow with the usual random- 
walk exponent 1/2, in agreement with experiments on two-dimensional arches pEj - 
Our observations on horizontal extensions of three-dimensional bridges have yielded 
a non-trivial exponent u^ n « 0.66. Putting all of this together, our results predict 
that long linear bridges are domelike; also, they are vertically diffusive but horizontally 
superdiffusive. Evidently, jamming in a three-dimensional hopper would be caused by 
the planar projection of such a dome. 

We now compare the results of this simple theory with data on bridge structures 
obtained from independent numerical simulations of shaken hard sphere packings |lUj . 
Figure El confirms that the mean angle is Gaussian to a good approximation, while 
Figure |H1 shows the measured size dependence of the variance (@ 2 )(s). The numerical 
data are found to agree well with a common fit to the first (stationary) term of (|2.1()jl 
- the 'transient' effects of the second term of ()2.10|) are too small to be significant at 
our present accuracy. We thus conclude that our simple theory captures the principal 
structural features of linear bridges. 



2.4- Discussion 

We end this section with the following remarks. First, more subtle effects, including the 
effects of transients via the second term of ()2.1Uj) . and the dependence of the parameters 
a 2 and a on the packing fraction $, are deserving of further investigation. Second, 



Competition and cooperation 

5 



10 



4 

^3 
® 

^2 
1 







• • o 


j. r\ rr n 

o $=(J.bb 
• $=0.58 






*\ 

8\ 
S\ 






. — 







0.2 



0.4 



0.6 



Figure 5. Plot of the normalised distribution of the mean angle (in radians) 
of linear bridges of size n = 4, for both volume fractions. The sinB Jacobian 
has been duly divided out, explaining thus the larger statistical errors at small 
angles. Full lines: common fit to (half) a Gaussian law. 
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Figure 6. Plot of the variance of the mean angle of a linear bridge, against 
size n, for both volume fractions. Full line: common fit to the first (stationary) 

0.093 and a = 0.55. The 'transient' effects of the 
are invisible with the present accuracy. 



term of (|2.10j) , yielding a^ q 
second term of (12.1 



we might expect that with increasing density $, branched structures would become 
more and more common; linear bridge formation, with its 'sequential' progressive 
attachment of independent blobs would then become more and more rare. Our theory 
should therefore cease to hold at a limit packing fraction $H m , which is qualitatively 
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reminiscent of the single-particle relaxation threshold density [S] (see section 3). Finally, 
our investigations suggest that long-lived bridges are natural indicators of sustained 
inhomogeneities in granular systems. 

3. On angles of repose: bistability and collapse 

Typically, the faces of a sandpile are inclined at a finite angle to the horizontal. This 
is the so-called 'angle of repose' in practice, it can take a range of values before 
spontaneous flow occurs, i.e., the sandpile becomes unstable to further deposition. The 
limiting value of this pre-avalanching angle is known as the maximal angle of stability 
9 m [Bj. Also, as a result of their athermal nature, sandpiles are strongly hysteretic; this 
results in bistability at the angle of repose [23 I2H], such that a sandpile can either 
be stable, or in motion, at any angle 9 such that #r < 9 < 9 m . Notwithstanding the 
above, it is possible for a sandpile to undergo spontaneous collapse to the horizontal; 
this is, in general, a rare event. We propose a theoretical explanation j7j below for both 
bistability at, and collapse through, the angle of repose via the coupling of fast and slow 
relaxational modes in a sandpile pQ. 

3.1. Coupled nonlinear equations: dilatancy vs. the angle of repose 

Our basic picture is that fluctuations of local density are the collective excitations 
responsible for stabilising the angle of repose, and for giving it its characteristic width, 

59 B = 9 m -9 R , (3.1) 

known as the Bagnold angle [221- Such density fluctuations may arise from, for instance, 
shape effects P or friction [SJ |2S] ; they are the manifestation in our model of Reynolds 
dilatancy |2Ej- 

The dynamics of the angle of repose 9{t) and of the density fluctuations (fi(t) are 
described jZj by the following stochastic equations, which couple their time derivatives 
9 and <j>: 

9 = - a 9 + b(f) 2 + AiTfc^), (3.2) 
0= -c</> + A 2 77 2 (t). (3.3) 

The parameters a, A 2 are phenomeno logical constants, while 7/i(t), T)2{t) are two 
independent white noises such that 

( Vl (t) Vj (t')) = 2 5 l] 5(t-t'). (3.4) 

The first terms in ()3.2j) and ()3.3j) suggest that neither the angle of repose nor 
the dilatancy is allowed to be arbitrarily large for a stable system. The second term 
in ()3.2j) affirms that dilatancy underlies the phenomenon of the angle of repose; in the 
absence of noise, density fluctuations constitute this angle. The term proportional to 
(f) 2 is written on symmetry grounds, since the the magnitude (rather than the sign) of 
density fluctuations should determine the width of the angle of repose. The noise in ()3.2j) 
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represents external vibration, while that in (|3.3|) embodies Edwards' compactivity [14j, 
being related to purely density- driven effects. We note that these equations bear more 
than a passing resemblance to those in the previous section on orientational statistics 
of bridges: the underlying reason for this similarity is the idea jSJ [7] that bridges form 
by initially aligning themselves at the angle of repose in a sandpile. 

Examining the above equations, we quickly distinguish two regimes. When the 
material is weakly dilatant (c 3> a), so that density fluctuations decay quickly to zero 
(and hence can be neglected), the angle of repose 9(t) relaxes exponentially fast to an 
equilibrium state, whose variance 

A 2 

2 eq = (3-5) 

is just the zero-dilatancy variance of 9. The opposite limit, where c -C a, and density 
fluctuations are long-lived, will be our regime of interest here. When, additionally, A x 
is small, the angle of repose has a slow dynamics reflective of the slowly evolving density 
fluctuations. These conditions can be written more precisely as 



The parameter 7, which sets the separation of the fast and slow timescales, is 

an inverse measure of dilatancy in the granular medium; small values of this imply 

a granular medium that is 'stiff' to deformation, resulting from the persistence of 

density fluctuations. The parameter e measures the ratio of fluctuations about the 

(zero-dilatancy) angle of repose to its full value in the presence of density fluctuations: 

from this we can already infer that it is a measure of the ratio of the external vibrations 

A 2 

to density- driven effects, which are explicitly contained in the ratio (~rj)- Realising that 



external vibrations and density /compactivity respectively drive fast and slow dynamical 
processes in a granular system, we see that a quantity which measures their ratio has all 
the characteristics of an effective temperature pQ in the slow dynamical regime of interest 
to us here. This temperature-like aspect will become much more vivid subsequently, 
when we discuss the issue of sandpile collapse. 

To recapitulate: the regime (j3.6|) that we will discuss below is characterised as 
low-temperature and strongly dilatant, governed as it is by the slow dynamics of density 
fluctuations. 

3.2. Bistability within 58^ : how dilatancy 'fattens 7 the angle of repose 




(3.6) 



(3.7) 




Suppose that a sandpile is created in regime (|3.6|) with very large initial values for the 
angle 9q and dilatancy 4>q. In the initial transient stages, the noises have negligible effect 
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and the decay is governed by the deterministic parts of (|3.2|) and (|3.3j) : 

9(t) = (9 -9 m )e- at + 9 m e- 2ct , (3.8) 
0(t) = o e- ct , (3.9) 

with 

«,„ « (3.10) 
a 

Thus, density fluctuations 4>(t) relax exponentially, while the trajectory 9(t) has two 
separate modes of relaxation. First, there is a fast (inertial) decay in 9(t) m 9q e~ at , until 
9{t) is of the order of # m ; this is followed by a slow (collective) decay in 9{t) ~ 9 m e~ 2ct . 
When (p{t) and 9{t) are small enough [i.e., <fi(t) ~ eq and 9{t) ~ #r, cf. (j3.11j) and (J3.13|) ] 
for the noises to have an appreciable effect, the above analysis is no longer valid. 
The system then reaches the equilibrium state of the full non-linear stochastic process 
represented by ()3.2|) and ()3.3|) . a full analytical solution of which is presented in 

In order to get a feeling for the more qualitative features of the equilibrium state, 
we note first that the equilibrium variance of <j)(t) is: 

< = — • (3-11) 

We see next that to a good approximation, the angle 9 adapts instantaneously to the 
dynamics of <fi(t) in regime (|3.6|) : 

9(t) « (3.12) 
a 

The two above statements together imply that the distribution of the angle 9{t) is 
approximately that of the square of a Gaussian variable. The typically observed angle 
of repose 6*r is the time-averaged value 

b<tL bA 2 , 

/eq 



9r = Wee = ^ = — * (3-13) 



ac 

Equation ()3.12|) then reads 

9(t)^9 R ^f. (3.14) 

Equation ()3.14j) entirely explains the physics behind the multivalued and history- 
dependent nature of the angle of repose [23 |2H]- Its instantaneous value depends 
directly on the instantaneous value of the dilatancy; its maximal (stable) value 9 m 
is noise-independent [cf. ()3.1()j) ] and depends only on the maximal value of dilatancy 
that a given material can sustain stably Sandpiles constructed above this will first 
decay quickly to it; they will then decay more slowly to a 'typical' angle of repose 
The ratio of these angles is given by 

r = #• (3 - 15) 

9r eq 

so that 9 m ^> 6*r for O ^> eq . Within the Bagnold angle 59b, (i- e - for sandpile 
inclinations which lie in the range #r < 9 < 9 m ), our simple theory also demonstrates 
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the presence of bistability. Thus, sandpiles submitted to low noise are stable in this 
range of angles (at least for long times ~ 1/c) ; on the other hand, sandpiles submitted 
to high noise (such that the effects of dilatancy become negligible in ()3.2j) ) continue to 
decay rapidly in this range of angles, becoming nearly horizontal at short times ~ 1/a. 

Our conclusions are that bistability at the angle of repose is a natural consequence 
of applied noise (tilt [2H] or vibration) in granular systems. For sandpile inclinations 9 
within the range 59b, sandpile history is all- important: depending on this, a sandpile 
can either be at rest or in motion at the same angle of repose. 

3.3. When sandpiles collapse: rare events, activated processes and the topology of 
rough landscapes 

When sandpiles are subjected to low noise for a sufficiently long time, they can collapse 
PQ, such that the angle 8(t) vanishes. Such an event is expected to be very rare in the 
regime (jH.fij) ; in fact it occurs only if the noise rji(t) in (J3.2)) is sufficiently negative for 
sufficiently long to compensate for the strictly positive term hcf 2 . It can be shown [Jj that 
the equilibrium probability for 9 to be negative, II = Prob(# < 0), scales throughout 
regime ()3.6|) as: 

n ^(2 £ ) 1 / 4 7 bAl 

n "iW4)^ c) ' c = 7^ = ^a 1 (3 ' 16) 

The scaling function .F(C) decays [Zj monotonically from J^(0) = 1 to JF(oo) = 0; to find 
out when the angle of repose first crosses zero, we should explore the latter limit, i.e. 
the regime ( > 1. Here, the equilibrium probability of collapse vanishes exponentially 
fast: 

n~expl~l^l I. (3.17) 



2 



The above suggests that sandpile collapse is an activated process, with a competition 
between 'temperature' e and 'barrier height' r y 2 . Collapse events occur at Poissonian 
times, with an exponentially large characteristic time given by an Arrhenius law: 

r~l/n~exp^^) j. (3.18) 

The stretched exponential with a fractional power of the usual 'barrier-height-to- 
temperature ratio' 7 2 /e is suggestive of glassy dynamics [3]; it also reinforces the idea 
that sandpile collapse is a rare event. 

While the reader is referred to a longer paper [7] for the derivation of the stretched 
exponential, the physics behind it is readily understood by means of an exact analogy 
with the problem of random trapping |29j . which we outline below. 

Consider a Brownian particle in one dimension, diffusing (with diffusion constant 
D) among a concentration c of Poissonian traps. Once a trap is reached, the particle 
ceases to exist, so that its survival probability S(t) is also the probability that it has not 
encountered a trap until time t. Assuming a uniform distribution of starting points, the 
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fall-off of this probability can be estimated by first computing the probability of finding 
a large region of length L without traps, and then weighting this with the probability 
that a Brownian particle survives within it for a long time t: 

r°° ( Tr 2 Dt\ 
S{t) ~ / exp -cL - dL. (3.19) 



The first exponential factor exp(— cL) is the probability that a region of length L is free 
of traps, whereas the second exponential factor is the asymptotic survival probability 
of a Brownian particle in such a region, exp(— Dq 2 t). The integral is dominated by a 
f27r 2 Dt\ 1/3 

saddle-point at L m \ , whence we recover the well-known estimate 

V c / 

3 ( n 2 



S(t) ~ exp \-- (2n 2 c 2 Dt) j . (3.20) 

Notice the similarity in the forms of ()3.17j) and ()3.2()|) : it turns out that the steps 
in their derivations are identical [7j, and form the basis of an exact analogy. In turn 
the analogy allows us to formulate an optimisation-based approach to sandpile collapse, 
which makes for a much more intuitive grasp of its physics. 

Accordingly, let us visualise the angle 9 as an 'exciton' whose 'energy levels' are 
determined by the magnitude of 9. It diffuses with temperature e in a frozen landscape 
of 4> (dilatancy) barriers of typical energy 7. Only if it succeeds in finding an unusually 
low barrier can it escape via ()3.17j) . to reach its ground state (9 = 0)- this of course 
corresponds to sandpile collapse. Taking the analogy a step further, we visualise the 
exciton as 'flying' at a 'height' 9, surrounded by 0-peaks of typical 'height' 7 in a rough 
landscape. Flying too low would cause the 9 exciton to hit a barrier fast, while flying 
too high would cause the exciton to miss the odd low barrier. It turns out [7j that flying 
at 9 ~ e 1 / 3 allows the exciton to escape via ()3.17j) (cf. the arguments leading to L ~ i 1 / 3 
above). Translating back to the scenario of sandpile angles, the above arguments imply 
the following: angles of repose that are too low are unsustainable for any length of 
time, given dilatancy effects, while angles that are too large will resist collapse. Thus 
optimal angles for sandpile collapse are found to scale as 9 ~ e 1 / 3 ; sandpiles with these 
inclinations show a finite, if small, tendency to collapse via ()3.17|) . 

Clearly, the frequency of collapse will depend on the topology of the (^-landscape; 
the form ()3.17j) was valid for a landscape with Gaussian roughness [Zj. What if the 
landscape is much rougher or smoother than this? - to answer this question, we look at 
two opposite extremes of non-Gaussianness. 

First, let us assume that density fluctuations are peaked around zero; typical 
barriers are low, and the ^-landscape is much flatter than Gaussian. The exciton's 
escape probability ought now to be greatly increased. This is in fact the case [Zj; it can 
be shown that in the 7 — > limit, the collapse probability scales as e 1 ^. Switching back 
to the language of sandpiles, this limit corresponds to a nearly non-dilatant material; 
it results in a 'liquid-like' scenario of frequent collapse, where a finite angle of repose is 
hard to sustain under any circumstances. 
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In the opposite limit of an extremely rough energy landscape, where large values 
of (p are more frequent than in the Gaussian distribution, one might expect the escape 
probability of the 9 exciton to be greatly reduced. If, for example, the jaggedness of 
the landscape is such that \4>(t)\ is always larger than some threshold the stretched 
exponential in ()3.17j) reverts (in the e C 1 regime considered) to an Arrhenius law in 
its usual form: 



In the language of sandpiles, this limit corresponds to strongly dilatant material; here, 
as one might expect, sandpile collapse is even more strongly inhibited than in (J3.17|) . 
Wet sand for example, is strongly dilatant; its angles of repose can be far steeper than 
usual, and still resist collapse. 

3.4- Discussion 

The essence of our theory above is that dilatancy is responsible for the existence of the 
angle of repose in a sandpile. We claim further that bistability at the angle of repose 
results from the difference between out-of-equilibrium and equilibrated dilatancies. We 
are also able to provide an analytical confirmation of the following everyday observation: 
weakly dilatant sandpiles collapse easily, while strongly dilatant ones bounce back. 

4. Compaction of disordered grains in the jamming limit: sand on random 
graphs 

Granular compaction is characterised by a competition between fast and slow degrees 
of freedom [T]; when voids abound, individual grains can quickly move into these to 
suit their convenience. As the jamming limit is approached, clusters of grains need to 
rearrange in order better to fill remaining partial voids: this process is necessarily slow, 
and eventually leads to dynamical arrest (3]. 

The modelling of granular compaction has been the subject of considerable effort. 
Early simulations of shaken hard sphere packings ^U], carried out in close symbiosis 
with experiment [SHI, were followed by lattice-based theoretical models jSU E2]; the 
latter could not, of course, incorporate the reality of a disordered substrate. Mean- 
field models [33] which could incorporate such disorder could not, on the other hand, 
impose the finite connectivity of grains included in Refs. [TO] I3T] 132*] . It was to answer 
the need of an analytically tractable model which incorporated finitely connected grains 
on fully disordered substrates, that we introduced jH] random graph models of granular 
compaction. 

Why random graphs? First, random graphs j33] are the simplest structures with 
a finite number of neighbours. Finite connectivity is a key property of grains in a 
sandpile with a gamut of consequences, ranging from kinetic constraints [33] to cascade 
processes in granular compaction |H]. Second, random graphs are among the simplest 




(3.21) 
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fully disordered constructs where, despite the existence of defined neighbourhoods of 
a site, no global symmetries exist. Disorder is an equally key feature of granular 
matter, even at the highest densities; its consequences include the presence of a 
range of coordination numbers jTU] for any sandpile, corresponding to locally varying 
neighbourhoods of individual grains, a feature which can be incorporated via locally 
fluctuating connectivities [H] in random graphs. 

Having motivated our choice of random graphs as a basis, we proceed below to 
describe our spin model of granular compaction. 

4-1. The three-spin model: frustration, metastability and slow dynamics 

The guiding factor in our choice of spin model is that it be the simplest model with 
frustration, metastability and slow dynamics; we will discuss the last two later, but 
remark at the outset that geometrical frustration is crucial to any study of granular 
matter. This concerns the fruitless competition between grains which try - and fail - to 
fill voids in the jamming limit, due either to geometric constraints on their mobility, or 
because of incompatibilities in shape or size. Our way of modelling this is via multi-spin 
interactions on plaquettes on a random graph [S]. We choose in particular a 3-spin 
Hamiltonian on a random graph (see Fig. E|) where N binary spins Si = ±1 interact in 
triplets: 

H = -pN = - C^SiSjSu (4.1) 

i<j<k 

Here, the variable C^k = 1 with i < j < k denotes the presence of a plaquette 
connecting sites i,j,k, while C^k = denotes its absence. Choosing C^k = 1(0) 
randomly with probability 2c/N 2 (1 — 2c/N 2 ) results in a random graph, where the 
number of plaquettes connected to a site is distributed with a Poisson distribution of 
average c - this exemplifies the locally varying connectivities mentioned above. The 
connection with granular compaction is made in accordance with Edwards' hypothesis 
[T3] : we interpret the local contribution to the energy in different configurations of the 
spins as the volume occupied by grains in different local orientations. 

This Hamiltonian has been studied on a random graph in various contexts jHEllSZI- 
It has a trivial ground state where all spins point up and all plaquettes are in the 
configuration + + + giving a contribution of —1 to the energy. Yet, locally, plaquettes 

of the type h, — I — , H (satisfied plaquettes) also give the same contribution, 

although one may not be able to cover the entire graph with these four types of plaquettes 
in equal proportions. This degeneracy of the four configurations of plaquettes with 
SiSjSk = 1 results in frustration, via the competitionbetweeia satisfying plaquettes locally 
and globally. In the former case, all states with even parity may be used, resulting in 
a large entropy and in the latter, only the + + + state may be used. Such frustration 
eventually leads to slow dynamics. 

This mechanism has a suggestive analogy in the concept of geometrical frustration of 
granular matter, if we think of plaquettes as granular clusters. When grains are shaken, 




Figure 7. A part of a random graph with triplets of sites forming plaquettes 
illustratng its local treelike structure (no planarity or geometric sense of 
distance are implied). 



they rearrange locally, but locally dense configurations can be mutually incompatible. 
Voids could appear between densely packed clusters due to mutually incompatible grain 
orientations between neighbouring clusters. The process of compaction in granular media 
consists of a competition between the compaction of local clusters and the minimisation 
of voids globally §. 

Another key feature of this model is the existence of metastable states. We note 
from Figure IHlthat two spin flips are required to take a given plaquette from one satisfied 
configuration to another; an energy barrier thus has to be crossed in any intermediate 
step between two satisfied configurations. This has a mirror image in the context of 
granular dynamics, where compaction follows a temporary dilation; for example, a 
grain could form an unstable ('loose') bridge with other grains before it collapses into an 
available void beneath the latter [21^1]]. This mechanism, by which an energy barrier has 
to be crossed in going from one metastable state to another, is an important ingredient 
in models of granular compaction [38J. 

4-2. How we tap the spins - dilation and quench phases 

Our tapping algorithm is a simplified version of the tapping dynamics used in cooperative 
Monte Carlo simulations of sphere shaking [ID] . We treat each tap as consisting of two 
phases. First, during the dilation phase, grains are provided with free volume to move 
into; next, in the quench phase, they are allowed to relax until a mechanically stable 

§ There are indications |H8I that global minimisation of voids can sometimes win out, when 
granular media are subjected to prolonged low-intensity vibration; this results in a first-order jump 
to a 'crystalline' state over significant length scales in the sample. 
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Figure 8. The phase space of three spins connected by a single plaquette. 
Configurations of energy — 1 (the plaquette is satisfied) are indicated by a black 
dot, those of energy +1 (the plaquette is unsatisfied) are indicated by a white 
dot 

configuration is reached. 

More technically, the dilation phase is modelled by a single sequential Monte Carlo 
sweep of the system at a dimensionless temperature T. A site i is chosen at random and 
nipped with probability 1 if its spin s, is antiparallel to its local field hi, with probability 
exp(— hi/T) if it is not, and with probability 0.5 if hi = 0. This procedure is repeated N 
times. Sites with a large absolute value of the local field hi thus have a low probability 
of flipping into the direction against the field; such spins may be thought of as being 
highly constrained by their neighbours. The dynamics of our 'thermal' dilation phase 
differs from the 'zero-temperature' dynamics used in jlU] where a certain fraction of 
spins is flipped regardless of the value of their local field. Our choice jH] reflects the 
following physics: if grains are densely packed ('strongly bonded' to their neighbours), 
they are unlikely to be displaced during the dilation phase of vibration. 

The grains are then allowed to relax via a T = quench, which lasts until the 
system has reached a blocked configuration [IT], where each site % has Sj = sgn(/ij) or 
hi = 0: thus, each grain is either aligned with its local field, or it is a 'rattler' [12] • Thus, 
at the end of each tap (dilation + quench), the system will be in a physically stable 
configuration [TU] . 
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4-3. The compaction curve 

Among the most important of our results is the compaction curve obtained by tapping 
our model granular medium for long times. This is shown in Figure El where three 
regimes of the dynamics can be identified. In the first regime, fast individual dynamics 
predominates, while in the second, one sees a logarithmic growth of the density via slow 
collective dynamics. The last regime consists of system-spanning density fluctuations in 
the jamming limit, where our quantitative agreement with experiment [IB] allows us to 
propose a cascade theory of compaction during jamming. 
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Figure 9. Compaction curve at connectivity c = 3 for a system of 10 4 spins 
(one spin is flipped at random per tap). The data stem from a single run with 
random initial conditions and the fit (dashed line) follows l|4.2|) with parameters 
Poo = .971, po = .840, D = 2.76, and r = 1510. The long-dashed line (top) 
indicates the approximate density 0.954 at which the dynamical transition 
occurs, the long-dashed line (bottom) indicates the approximate density 0.835 
at which the fast dynamics stops, the single-particle relaxation threshold. 



4-3.1. Fast dynamics till SPRT: each grain for itself! At the end of the first tap, 
each grain is connected to more unfrustrated than frustrated clusters. This is a direct 
result of the first tap being a zero-temperature quench: any site where this was not the 
case would simply flip its spin. More generally, a fast dynamics occurs in this regime 
whereby single grains locally adopt the orientation that, finally, optimises their density; 
this density po has been termed [Hj the single-particle relaxation threshold (SPRT). 

If we neglect correlations between the local fields of neighbouring sites, we can 
arrive at a simple population dynamics model of the fast dynamics, details of which 
can be found in [Sj. It yields a value of the SPRT, p = 0.835 (shown as a dotted line 
in Figure El) which is much higher than the value of the density p(= 0.49) of a typical 
blocked configuration. 




Competition and cooperation 



21 



This is quite an extraordinary result; it implies that despite the exponential 
dominance of blocked configurations, random initial conditions preferentially select a 
higher density corresponding to the SPRT pp. This prediction of an overshoot in the 
density achieved by fast dynamics has also, strikingly, been confirmed in independent 
lattice-based models of granular compaction, and points to a strong non-ergodicity 
in the fast dynamics of individual grains. We will discuss this further in Section 5.3. 

Another significant feature of this regime is that a fraction of spins is left with 
local fields exactly equal to zero, which thus keep changing orientation |44.. These 
are manifestations in our model of 'rattlers' [12], i.e. grains which keep changing their 
orientation within well-defined clusters jTHj. They will later (cf. Section 5.3) be used as 
a tool to probe the statistics of blocked configurations jH] . 

To summarize: each grain reaches its locally optimal configuration via fast 
individual dynamics, resulting in the attainment of the SPRT density. All dynamics 
after this point is perforce collective. 

4-3.2. Slow dynamics of granular clusters: logarithmic compaction. The second stage 
of our tapping dynamics is fully collective: it removes some of the remaining frustrated 
plaquettes as clusters slowly rearrange themselves. A logarithmically slow compaction 
results [HOI EH], leading from the SPRT density po to the asymptotic density p^. The 
resulting compaction curve may be fitted, with D and r being characteristic constants, 
to the well-known logarithmic law [30J: 

pit) = poo - (poo - p )/(l + 1/D ln(l + t/r)) , (4.2) 

This can be written more transparently as 1 + t(p)/r = exp{D— — — } , a form 

Poo - P 

which makes clear that the dynamics becomes slow (logarithmic) as soon as the density 
reaches po- Although most grains are firmly held in place by their neighbours in this 
regime, cascade-like changes of orientation can occur. For example, if some grains 
change orientation during the dilation phase, this would change the constraints on their 
neighbours; importantly, the freer dynamics of rattlers could also alter local fields in 
their neighbourhood, and cause previously blocked grains to reorient. Reorientation 
in cascades jH] would then ensue, leading to collective granular compaction upto the 
asymptotic density p^. We identify this with the density of random close packing [IS] 
(RCP) and associate it with a dynamical phase transition (HUE]. 

4-3.3. Cascades at the dynamical transition. With increasing density, free-energy 
barriers rise up causing the dynamics to slow down according to ()4.2|) . The point where 
the height of these barriers scales with the system size marks a breaking of the ergodicity 
of the dynamics; an exponential number of valleys appears in the free-energy landscape. 
Our theoretical |Hj value for the dynamical transition, shown as a horizontal line in 
Figure 01 is in good agreement with the asymptotic density p^ reached (numerically) 
by the tapping dynamics. 
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Also shown in Figure El are marked fluctuations around the logarithmic compaction 
law, especially as the jamming limit is approached; their correlations over several octaves 
have been the subject of detailed experimental investigations 43J. To compare the 
results of our model with experiment, we follow [43 and use Fourier transforms to plot 
the power spectrum of the timeseries p(t) (in different frequency bands) against time in 
Figure [TUl 
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Figure 10. The density fluctuations as a function of time resulting from 
1024 taps are plotted as the topmost trace. The successive plots are of the 
power spectrum against time, in different frequency octaves. The power in the 
first octave (frequency 1/(1024 taps) -2/(1024 taps)) is the bottommost trace, 
second octave (frequency 2/(1024 taps) - 4/(1024 taps)), above it, and so on 
to the top. Note that the fluctuations of the power in the different frequency 
bands are strongly correlated; they correspond to sudden changes in the density 
(topmost trace). 

Our results indicate, as in the experimental data, that there are 'bursts' in the 
power spectrum fluctuations: our decomposition of these bursts over several octaves 
shows that they are caused by strong correlations of noise power over a wide range 
of frequencies. Importantly, the correlation matrices we obtain || are in quantitative 
agreement with experiment [4*3*] . 

In our model, we can trace such bursts as being due to cascades of spin-flips. As 
mentioned before, these arise from the change in local fields caused by the flipping 
of a single spin (or several spins); this instability propagates through ever larger 
neighbourhoods, causing correlated bursts in noise power fluctuations. Note that this 
cascade mechanism is entirely absent from generic fully connected models |33], where 
each spin interacts with all spins in the system. Interestingly, it is also absent from 
the finitely-connected parking- lot model [48]; this is because the creation and filling of 
'parking lots' by 'cars' in that model does not cause the appearance of further parking 
lots, in the way spins flips may trigger a cascade. Put another way, the absence of 
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competition between individual and collective interactions in the parking-lot model 
results in an absence of cascades there, despite its finite connectivity. 

We use the above insights to argue [Hj that in real granular media, the observed 
correlations [HJ| of density fluctuations are due to a cascade process in granular 
compaction near the jamming limit. Here, orient at ional /positional changes in strongly 
constrained grains give rise to propagating instabilities, leading to a near-global 
rearrangement of the granular medium. Pictorially, the movement of a single grain 
in this regime is only possible as the consequence of a system-wide cooperative motion 
of grains: this leads to sharp changes in overall density, and to the observed 'bursts' in 
the power spectrum of density fluctuations 43J. 

4-4- Realistic amplitude cycling; how granular media jam at densities lower than 
close-packed. 

Our random graphs model has also been used to simulate amplitude cycling, an 
experimental [30] protocol on tapped granular media. Here, the granular medium is 
tapped at a given amplitude T for a time r, after which its amplitude is changed by an 
infinitesimal 5T; this process is repeated for cycles of increasing and decreasing T. The 
control parameter turns out to be the so called 'ramp rate', which is the ratio ST/t; 
this is a measure [31] of the 'equilibration' allowed to the granular medium. Clearly low 
values of this will correspond to quasi-static processes, while large values will correspond 
to near-adiabaticity. 

The most simple-minded application of this protocol in our model results in the 
scenario of Figure El At both high and low cycling rates, the density p first reaches 
the SPRT p , increases with increasing amplitude, and decreases again at large values 
of T. Thereafter, p always decreases with increasing T. The part of the curve where T 
is increased for the first time has been termed [3D] the irreversible branch; the reversible 
branch refers to the trajectory traced out by all successive increases and decreases of 
tapping amplitude. The results of Figure ^3 i n agreement with many other models 
[3U I3E1 HH] j suggest that as the ramp rate is decreased, the system will eventually attain 
the RCP density p^. In particular, these models predict that in the limit of near-zero 
ramp rates, the irreversible branch disappears, with p becoming a single- valued function 

of r. 

This prediction is in direct contradiction to the experimental results of [30]; these 
suggest that at least for experimentally realisable times, low-amplitude shaking does not 
result in RCP being reached. Instead, they suggest that some grains in the jamming 
limit of a granular assembly are so strongly constrained that they will never be displaced 
by low-amplitude taps. 

In order to model the above experimental scenario, we modify the simplest picture 
of amplitude cycling presented above. First, we realise that in the dynamics of our model 
described thus far, sites with a high local field (corresponding to strongly constrained 
grains) may be flipped at any finite value of V with a correspondingly small but finite 
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Figure 11. Amplitude cycling: T is varied in both directions between to 2 at 
two rates <5r = 10~ 4 , 10 -5 per tap. with r = 1. The lower ramp rate (shown 
by the dotted line) results in a higher final density. 

probability. This is what eventually leads the system to the RCP density p^. 

To prevent this drift to RCP, one could naively think of introducing a threshold 
in local fields such that spins with fields above this threshold are not nipped. It turns 
out |H], however, that this is insufficient; the orientational dynamics of neighbouring 
spins will always loosen the constraints on previously blocked spins in the end, and lower 
their local fields below any given threshold. The above implies that the constraints 
on grains are not related to orientational frustration alone; it was suggested jH] that 
they might also be mechanical in nature, related to force networks [Tj\ between grains. 
Further, it seemed reasonable to suppose that such effectively immobile [HI] 'blockages' 
could only be removed kinetically by the imposition of large intensity (T) vibrations. 

Accordingly, the concept of low-amplitude pinning of grains was introduced |Hj: 
assign to each site i a real number between zero and one, such that only grains with 
?"j < T (mechanically constraining forces less than external vibration intensity) will be 
free to move. This modification could in principle lead to a lower value of p^ after 
amplitude recycling: for example, spin plaquettes (granular clusters) generated during 
the high-amplitude part of the cycle, would be effectively immobile at lower amplitudes, 
leading to wasted space. 

However, this too is insufficient to stop the evolution of the system to poo. It turns 
out 8] that jamming at lower densities can only be achieved if low-amplitude pinning is 
combined with the choice of an extremely low ramp rate, via a large 'equilibration time' 
t. If this is done, grains are allowed to 'equilibrate' at each amplitude, thus making 
sure that a steady-state of the density is always reached. This leads, finally, to the 
low-amplitude immobilisation of granular clusters created at high amplitudes, to the 
consequent blocking of voids, and hence to 'jamming' [3UJ ED] at densities lower than 
RCP (see Figure EJ) 
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Figure 12. The asymptotic density for tapping amplitudes ranging from 
T = 0.2 to T = 2 in steps of 0.2. The density measured after 10 7 taps at 
each amplitude and convergence to a steady-state each time was checked. 



Our results jHj on modified amplitude cycling thus indicate that it is important to 
have mechanical pinning as well as long equilibration times j22] for jamming to occur. 
In fact, the random configurations of immobile spins at each value of T can be viewed jH] 
as an additional quenched disorder, and their effect on neighbouring mobile spins, as 
an additional random local field. Our results demonstrate also a rather fundamental 
difference between excitations in glassy systems and granular media. In glasses, one 
would expect the configurations of spins reached at high values of temperature and 
subsequently frozen, immediately to alter the behaviour of the system at lower values 
of temperature. In athermal media like granular systems, however, it is important to 
allow equilibration at each value of the shaking intensity T, in order even to begin to 
observe the hysteretic effects which result in jamming at lower-than-RCP densities. 

4-5. Discussion 

In the above, we have presented a finitely connected spin model on random graphs jH], 
with the aim of examining the compaction of tapped granular media. Fast non-ergodic 
relaxation of individual grains terminated at the SPRT density, in our compaction 
curve; collective relaxation followed, manifest first by logarithmic compaction, and next 
by system-wide density fluctuations around RCP, both of which match experimental 
results |Sni HH]- Our model explains the latter in terms of a cascade process that 
occurs near the jamming limit of granular matter. Also, in contradistinction to other 
models jHU EHl HH], our results indicate that jamming at densities lower than RCP 
occurs as a result of competition between mechanical and orientational frustration, during 
amplitude cycling. We leave the discussion of the configurational entropies jHl 153] 
generated by our tapping algorithm to the concluding section of this article. 



Competition and cooperation 



26 



5. Shape matters 

In the concluding section of this review, we explore the effect of grain shapes in 
granular compaction. Our model is based on the following picture. Consider a box 
of sand in the presence of gravity; for ease of visualisation, we think of this as being 
constituted of vertical columns. In the jamming limit, diffusion of grains between 
columns is inhibited [3T] due to the absence of holes (grain-sized voids); compaction 
proceeds instead by grain reorientation within each column |Hj to minimise the size of 
the partial voids that persist. This intuitive picture is verified by results of computer 
simulations ^Hj which show that correlations in the transverse (inter-column) direction 
are negligible compared to those in the longitudinal (intra-column) direction. We thus 
focus on a column model of grains in the jamming limit 0. 

Each ordered grain occupies one unit of space, while each disordered grain occupies 
l+e units of space, with e a measure of the partial void trapped by mis orientation. A 
reorientation of the grain to an ordered ('space-saving') state frees up the partial void, 
for use by other grains to reorient themselves. This cascade-like picture of compaction in 
the jamming limit resembles that presented in the previous section, where random graph 
models [Hj were discussed. Also, as there, the response of grains to external dynamics 
is via the local minimisation of void space modelled by a local field. 

5.1. A model of regular and irregular grains 

In our column model, grains are indexed by their depth n measured from the free surface. 
Each grain can be in one of two orientational states - ordered (+) or disordered (— ) 
- the 'spin' variables {a n = ±1} thus uniquely defining a configuration. Exactly as in 
the random graphs model |HJ presented above, a local field h n constrains the temporal 
evolution of spin a n , such that excess void space is minimised. 

In the presence of a vibration intensity T, grains reorient with an ease that depends 
on their depth n within the column (grains at the free surface must clearly be the 
freest to move!), as well as on the local void space h n available to them. The transition 
probabilities governing this are: 



The dynamical length £ dyn Ell effectively defines the boundary layer of the column; 
within this dynamics are fast, while well beyond it, they are slow. The local field h n is 
a measure of excess void space jHJ: 



where m+ and m~ are respectively the numbers of + and — grains above grain n. The 
definition Equation (|5.2|) is such that a transition from an ordered to a disordered state 
for grain n is hindered by the number of voids that are already above it, as might be 
expected for an ordering field in the jamming limit. 

In the r — > limit of zero-temperature dynamics [S], the probabilistic rules (JO]) 
become deterministic: the expression a n = sign h n (provided h n ^ 0) determines the 



w(a n = ± -> a n = =f) = exp(-n/£ dyn =f K/T). 




h n = em„ -m. 



(5.2) 
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ground states of the system. Frustration j2] manifests itself for e > 0, which leads to a 
rich ground-state structure, whose precise nature depends on whether e is rational or 
irrational. We mention for completeness that the case e < discussed in earlier work 
pTT] corresponds to a complete absence of frustration and a single ground state of ordered 
grains. 

For irrational e, no local field h n can ever be zero (cf. (|5.2Jl ). Noting that irrational 
values of e denote shape irregularity, we conclude that the excess void space is nonzero 
even in the ground state of jagged grains. Their ground state, far from being perfectly 
packed, turns out 9J to be quasiperiodic. 

Regularly shaped grains correspond to rational e = p/q, with p and q mutual primes. 
We see from (j5.2j) that now, some of the h n can vanish; these correspond, as noted in the 
previous section, to rattlers. A rattler at depth n thus has a perfectly packed column 
above it, so that it is free to choose its orientation [HI E3 E2j For regular grains in their 
ground state, rattlers occur periodically (as in crystalline packings!) at points such that 
n is a multiple of the period p + q.\\. Every ground state is thus a random sequence of 
two patterns of length p + q, each containing p ordered and q disordered grains; this 
degeneracy leads to a zero-temperature configurational entropy or ground-state entropy 
S = In 2/(p + q) per grain. 

5.2. Zero-temperature dynamics: (ir)retrievability of ground states, density 
fluctuations and anticorrelations 

Regular and irregular grains behave rather differently when submitted to zero- 
temperature dynamics. The (imperfect) but unique ground state for irregular grains 
is rapidly retrieved; the perfect (and degenerate) ground states for regular grains never 
are, resulting in density fluctuations. 

We recall the rule for zero-temperature dynamics: 

o~n > sign h n . (5.3) 

Starting with irregular grains (with a given irrational value of e) in an initially 
disordered state, one quickly recovers the ground state with zero-temperature dynamics. 
The ground state in fact propagates ballistically from the free surface to a depth 
L(t) fa V(e) t [0] at time t, while the rest of the system remains in its disordered initial 
state. When L(t) becomes comparable with £ dyn , the effects of the free surface begin 
to be damped. In particular for t ^> C,dyn/V(e) we recover the logarithmic coarsening 
law L(t) ~ £dynl n ^ also seen in other theoretical models [HI EH] of the slow relaxation 
of tapped granular media 30J. To recapitulate, the ground state for irregular grains 
is quickly (ballistically) recovered with zero-temperature dynamics, until the boundary 

|| For example, when e = 1/2, each disordered grain 'carries' a void half its size; units of perfect packing 

must be permutations of the triad H , where two 'half voids from each of disordered grains are 

perfectly filled by an ordered grain. The stepwise compacting dynamics [5| selects only two of these 
patterns, H and — I — 
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layer £d yn is reached; below this, the column is essentially frozen, and coarsens only 
logarithmically. 

For regular grains with rational e, the local field h n in (J5.3)) vanishes for rattlers. 
Their dynamics is stochastic even at zero temperature, since they have a choice of 
orientations: a simple way to update them is according to the rule a n — > ±1 with 
probability 1/2. This stochasticity results in an intriguing dynamics even well within 
the boundary layer £ dyn , while the dynamics for n ^> £d y n is as before logarithmically 
slow 

In what follows, we will focus on the fast dynamics within the boundary layer. Our 
main result is that zero-temperature dynamics does not drive the system to any of its 
degenerate ground states, but instead engenders a fast relaxation to a non-trivial steady 
state, independent of initial conditions, which consists of unbounded density fluctuations. 
This recalls density fluctuations close to the jamming limit [HUSH!, in other studies of 
granular compaction. 

5 
4 

CV! a 3 

£ 2 
1 



2 4 6 8 

In n 

Figure 13. Log-log plot of W% = (h^) against depth n, for zero-temperature 
dynamics with e = 1. Full line: numerical data. Dashed line: fit to asymptotic 
behaviour leading to Q5.4JI (after j^j). 




Fig. ^1 shows the variation of these density fluctuations as a function of depth n: 

Wl = {h 2 n ) « An 2/3 , A ps 0.83. (5.4) 

The fluctuations are approximately Gaussian, with a definite excess at small values: 
\h n \ ~ 1 W n . We recall that non-Gaussianness was also observed in experiments on 
density fluctuations in tapped granular media |Hi]; in our theory here, we interpret it 
in terms of grain (anti) correlations. If grain orientations were fully uncorrelated, one 
would have the simple result (h%) = ne, while ()5.4|) implies that (h%) grows much more 
slowly than n. 
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Figure 14. Scaling plot of the orientation correlation function c m ^ n for n ^ m 
in the zero-temperature steady state with e = 1, demonstrating the validity 
of (|5.5|) and showing a plot of (minus) the scaling function F (after j^j). 



It turns out that at least within a dynamical cluster of radius n 2 ^ 3 ) jU], the 

orientational displacements of each grain are fully anticorrelated. Fig. shows that 
the orientation correlations c m ^ n = (a m a n ) scale as [§] 

~A 1 T? { n-m \ . . 



w m w n \w m w n 

where the function F is such that f*™ F(x)dx = 1. We find also that, within such 
a dynamical cluster, the fluctuations of the orientational displacements are totally 
screened: Z) n ^m c m,n ~ — c m , m = —1. These results recall the anticorrelations in 
grain displacements observed in independent simulations of shaken hard spheres close 
to jamming [TQj ; there they corresponded to compaction via bridge collapse, as upper 
and lower grains in bridges [6J collapsed onto each other, releasing void space. 

5.3. Rugged entropic landscapes: Edwards' or not? 

The most remarkable feature of our column model is, arguably, the rugged landscape of 
microscopic configurations visited during the steady state of zero-temperature dynamics 
(for regular grains); this is all the more striking because the macroscopic entropy is flat, 
in agreement with Edwards' hypothesis |T4"| . 

The entropy of the steady state of zero-temperature dynamics is denned by the 
usual Boltzmann formula 

S=-5>(C)lnp(C), (5.6) 

c 

where p(C) is the probability that the system is in the orientation configuration C in the 
steady state, and the sum runs over all the 2 n configurations of a system of n grains. 
This can be estimated theoretically by using (|5.4j) . Consider n as a fictitious discrete 
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Figure 15. Plot of the measured entropy reduction AS in the zero-temperature 
steady state with e = 1, against n < 19. Symbols: numerical data, for t ~ 10 9 
and n « 20. Full line: fit AS = (62 Inn + 53)10~ 3 n 1 / 3 . 

time, with the local field h n as the position of a random walker at time n. For a free 
lattice random walk of n steps, one has (h%) = n; as all configurations are equiprobable. 
the entropy reads Sflat = nln2. For a column of regularly shaped grains, our model 
predicts instead (h%) = W% -C n; the entropy S of our random walker is therefore 
reduced with respect to Sfl at . The entropy reduction [S3j AS = 5fl a t — S 1 = nln2 — S 
can be estimated [9\ to be 

n I 

^ ~ E 7772- - - 1/3 - ( 5 - 7 ) 

m=l m 

Evaluating the steady-state entropy S numerically, using ()5.6|) and measuring all 
configurational probabilities p(C), we find (cf. Figure ITo^l that AS is small; for example 
for n = 12, we have AS ~ 0.479, in good agreement with the results of Eq. (|5.7|) . This 
is a convincing demonstration that anticorrelations (see previous section) lead to only 
microscopic corrections to the overall dynamical entropy of the steady state, which is 
flat, in agreement with Edwards' hypothesis [T4] . 

To investigate the effect of the constraints, we plot the normalised configurational 
probabilities 2 12 p(C) for a column of 12 grains against the 2 12 = 4096 configurations 
C in Figure HB1 Note that the actual values of the configurational probabilities p{C) 
are microscopically small! At this microscopic scale, however, the entropic landscape is 
startlingly rugged; some configurations are clearly visited far more often than others. It 
turns out that the most visited configurations are the ground states of the system (empty 
circles). We suggest that this behaviour is generic: i.e., the dynamics of compaction in 
the jammed state leads to a microscopic sampling of configuration space which is highly 
non-uniform, so that its ground states are visited most frequently. Our model thus 
provides a natural reconciliation between, on the one hand, the intuitive perception 
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Figure 16. Plot of the normalised probabilities 2 12 p(C) of the configurations 
of a column of 12 grains in the zero-temperature steady state with e = 1, 
against the configurations C in lexicographical order. The empty circles mark 
the 2 6 = 64 ground-state configurations, which turn out to be the most probable 
(after 0). 



that not all configurations can be equally visited during compaction in the jamming 
limit, that the most compact configurations should be the most visited; and, on the 
other, the flatness hypothesis of Edwards, which states that for large enough systems, 
the entropic landscape of visited configurations is flat [T4*] . 

The dynamical entropy generated by the random graphs model [Hj of the previous 
section is also reconcilable with Edwards' flatness ^3], at least in the jamming limit 
discussed above. This was explored via rattlers ( sites % such that the local field hi = 0) 
in the blocked configurations generated after each tap. We have seen above that they 
have a rather crucial role to play the density fluctuations of our column model [§]; it turns 
out that they are also a good probe of Edwards' flatness under the tapping dynamics of 
our random graphs model. If blocked states at a given density are equiprobable, a plot 
of the fraction of connected rattlers versus the density should reproduce this jHj- 

Figure El shows the results for four single runs of plotting the fraction of rattlers 
g against density p, at increasing amplitudes of vibration T. The dashed line and full 
lines correspond respectively to quenched and annealed replica symmetric averages for 
g, assuming Edwards' flatness [Sj. We notice that there is a reasonable congruence of 
all the numerical results and the (theoretically more accurate) quenched average at the 
asymptotic density poo. Thereafter, there are systematic divergences with lower density 
and higher T. 

We can draw the following conclusions from this. First, at the jamming limit 
near RCP, the dynamically generated entropies are flat, in accord with Edwards' 
hypothesis [Hj, as well as with the results of our column model Second, as we move 
to the regimes of higher vibration and lower density, the entropic landscape gets rougher 
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Figure 17. The fraction g of connected rattlers during 4 runs of a tapped 
random graph 8 with N = 1000, c = 3 at T = 0.4 (dots), T = 0.56 (+), T = 0.7 
(x), and T = 1.5 (circles). The solid and dashed lines correspond respectively 
to annealed and quenched theoretical values corresponding to Edwards' flatness. 
The vertical lines indicate the approximate values for po (left line) and /Jqo (right 
line) respectively. 

- one can imagine a process whereby the roughening visible on microscopic scales near 
jamming (cf. Figure ITT)|) begins to increase to macroscopic scales as one moves away from 
jamming. In this regime, we observe that configurations which are dynamically accessed 
by tapping (cf. the symbols in Figure 1X7)1 correspond to higher than typical densities 
(dashed and full lines in Figure 1X7)) - we recall from section 4.3.1 that this occurs when 
non-ergodic fast dynamics dominates configurational access. Putting all of this together, 
we conclude that also according to our random graphs model, Edwards' flatness [Xlj 
governs dynamically generated configurational entropies when slow (ergodic) dynamics 
predominate in the jamming limit; there are however, systematic deviations from flatness 
when fast dynamics predominate, in the regime of higher tapping amplitudes and lower 
densities. 

Configurational entropies of strongly nonequilibrium models with slow dynamics, 
are, however, not generically flat. To demonstrate this, we present results for a model 
of nonequilibrium aggregation, which despite its origins in cosmology (55] turns out to 
have applications in the gelation of stirred colloidal solutions [SHJ EZj ■ This 'winner- 
takes-all' model of cluster growth, whereby the largest cluster always wins, manifests 
both fast and slow dynamics. In mean field, the slow dynamical phase results in at 
most one surviving cluster at asymptotic times; however, on finite lattices, there can 
be many metastable clusters which survive forever, provided they are each isolated from 
the others (Figure IX5j) . 

We remark that this 'isolation' of surviving sites implies a very strong 
anticorrelation between neighbouring sites in this model; that is, each survivor must 
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I'igurc 18. A typical pattern of surviving clusters on the square lattice for 
the cluster aggregation model of Refs. jSHlEZI- Black (resp. white) squares 
represent a n = 1 (resp. a n = 0), i.e., surviving (resp. dead) sites. The left 
panel shows a 150 2 sample, while the right panel is enlarged (40 2 ) for clarity. 



have voids around it, or run the risk of dying out. These anticorrelations are manifest 
in Figure both for cluster survival and cluster mass on a one- dimensional version 
of the model. The presence of such anticorrelations and of competition between slow 
and fast dynamics in a nonequilibrium context, suggests strong analogies between this 
model joTH loTj and our random graphs [S] and column models. We might therefore 
naively expect some version of Edwards' flatness to hold; however, our results [HE] 
suggest that it does not. 

In conclusion, we emphasise that Edwards' flatness in the landscape of 
configurational entropies is not the generic fate of strongly nonequilibrium models with 
slow dynamics, even when they have many features in common. The similarity between 
our column model [H] and the random graphs model jH] discussed in the previous section 
is thus all the more remarkable; both models manifest Edwards' flatness in the jamming 
limit, deviating from it whenever free volume constraints are relaxed. 

5-4- Low-temperature dynamics along the column: intermittency 

We now return to the investigation of the column model pj, to do with its low- 
temperature dynamics. For rational e, the presence of a finite but low shaking 
intensity merely increases the magnitude of density fluctuations jSO], given that the 
zero-temperature dynamics is in any case stochastic. However, for irrational e, low- 
temperature dynamics introduces an intermittency in the position of a surface layer, 
this has recently been observed in experiments on vibrated granular beds 

This happens as follows: when the shaking amplitude T is such that it does not 
distinguish between a very small void h n and the strict absence of one, the site n 'looks 
like' a point of perfect packing. The grain at depth n then has the freedom to point the 
'wrong' way; we call such sites excitations, using the thermal analogy. The probability 




Figure 19. Plot of correlation functions for the cluster aggregation model 
of Refs. |56l I57j against the distance n along the chain. Empty symbols: 
correlation C a {n) of the survival index. Full symbols: correlation C x (n) of 
the reduced mass. 



of observing an excitation at site n scales as U(n) « exp(— 2\h n \/T). The uppermost site 
n such that \h n \ ~ T <C 1 will be the 'preferred' excitation; it is propagated ballistically 
(cf. zero-temperature irrational e dynamics) until another excitation is nucleated above 
it. Its instantaneous position Af(t) denotes the layer at which shape effects are lost in 
thermal noise, i.e., it separates an upper region of quasiperiodic ordering from a lower 
region of density fluctuations (j5.4|) . 

Fig. ED shows a typical sawtooth plot of the instantaneous depth J\f(t), for a 
temperature T = 0.003. The ordering length, defined as (A/"), is expected to diverge 
at low temperature, as excitations become more and more rare; we find in fact jU] a 
divergence of the ordering length at low temperature of the form (A/") ~ 1/ (T| In r|). 
This length is a kind of finite-temperature equivalent of the 'zero-temperature' length 
£dyn, as it divides an ordered boundary layer from a lower (bulk) disordered region. We 
emphasise once again that fast dynamics predominates in both zero-temperature and 
finite temperature boundary layers (£d y n and (A/")), with slow dynamics setting in for 
column depths beyond these. 

5.5. Discussion 

We have discussed the effect of shape in granular compaction near the jamming limit, 
via a column model of grains i9J. Our main conclusions are that jagged (irregular) 
grains are characterised by optimal ground states, which are easily retrievable, while 
smooth (regular) grains cannot retrieve their ground states of perfect packing; in the 
latter case, even zero-temperature dynamics results in density fluctuations. We predict 
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Figure 20. Plot of the instantaneous depth Af(t) of the ordered layer, for e = $ 
(the golden mean) and T = 0.003. Dashed lines: leading nucleation sites given 
by Fibonacci numbers (bottom to top: F\\ = 89, F\2 = 144, F13 = 233) 
(after 0). 



also that grain irregularities result in a surface intermittency for low- amplitude shaking, 
in agreement with recent experimental observations |5H] . 

6. Conclusions 

We have here reviewed some of our approaches to granular dynamics, now well known 
to consist of both fast and slow relaxational processes jl]. In the first case, grains 
typically compete with each other, while in the second, they cooperate. The dynamics 
of bridge formation is a typical result of cooperation, as is the relaxation of the angle of 
repose; competition between density fluctuations and external driving forces can, on the 
other hand, result in sandpile collapse. In the random graphs model presented above, 
the SPRT density separates regions of cooperation and competition, each with its own 
distinctive features. Finally, while slow dynamics predominate deep inside our column 
model of compacting grains, fast dynamics gives rise to strikingly rough configurational 
landscapes and surface intermittency. 
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